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Abstract 

I extend to QCD an efficient method for lattice gauge tlieory with dynamical 

^ , fermions. Once the eigenvalues of the Dirac operator and the density of states of pure 

gluonic configurations at a set of plaquette energies (proportional to the gauge action) 

C^ ' are computed, thermodynamical quantities deriving from the partition function can 

^^ , be obtained for arbitrary flavor number, quark masses and wide range of coupling 

I \ constants, without additional computational cost. Results for the chiral condensate 

^^ ' and gauge action are presented on the 10^ lattice at flavor number A'^^ = 0, 1, 2, 3, 

4 and many quark masses and coupling constants. New results in the chiral limit for 

the gauge action and its correlation with the chiral condensate, which are useful for 

p ,' analyzing the QCD chiral phase structure, are also provided. 



1 Introduction 



H \ Although it is believed to be the most powerful non-perturbative approach to quantum 
field theory, lattice Monte Carlo techniques still suffer from systematic errors. The sources 
of systematic error are mainly finite size effects, finite lattice spacing, quenched approxi- 
mation, and the use of unphysical quark masses. Finite volume effects and lattice spacing 
errors are relevant, and can be significantly reduced through the implementation of the 
Symanzik improvement [|l], |^, H, §, ||, |, |^, H- Quenched approximation ignores the feed- 
back of fermions (the determinant of the Dirac operator A in the measure is constraint to 



be det A = 1). More recent investigations @, |T0|, |Tl|] suggest that quenched approximation 



may lead to a systematic deviation from experiment. 

The most popular algorithms for generating unquenched configurations are the hybrid 
Monte Carlo (HMC) |]12| and/or hybrid molecular dynamics (HMD) algorithms JlSl. In 



these algorithms, each trial configuration is generated by solving the equation of motion 
with some number of steps, and in each step, numerical inversion of the Dirac matrix 



has to be done. For a different set of bare parameters of "sea" quark mass rrisea, flavor 
number Nf, and coupling constant g, new independent simulations are required. Tlierefore, 
tliese algoritlims are still very expensive and require liigfi performance supercomputers, in 
particular for small quark masses, even on a moderate lattice. 

Without loss of generality, let us discuss the unimproved lattice QCD with Kogut- 
Susskind quarks, which has the action S = Sg + Sf. Here 

S, = -|-EReTr(t/p), 

x,y 
Up = Uf,{x)U^{x + iJ.)Ul{x + u)Ul{x), 

1 r 

+ -77^ (x) \U^{x)5^^y_^ -Ulix- /U)(5x,y+^ 
T^^{X) = (_l)^i+^2+-+^M-l. 

(1) 

The bare coupling constant g is related to (3 hj (3 = 2Nc/g'^ with Nc = 3 the number of 
colors. The corresponding partition function is 

Z = f[dU] e-^« (det A(m, ^7))"^^/' . (2) 

In most simulations, people (e.g. [0) usually distinguish "sea" quarks (for the configu- 
rations) from the "valence" quarks (for the quark propagators), due to the heavy costs in 
generating the full QCD configurations at different set of sea quark masses. Let us take 
the chiral condensate as an example, which is meant to compute 

, T , ^ J[dU]TiA-\m,ai, U)e-'^ (det A(m,ea, U)f''^ 
VN, ![dU]e-^^ (det A(m,ea, U)f^/' 

(3) 

at a set of valence quark mass irtyau while always using a full configuration at a given sea 
quark mass irtsea 7^ ^T^vai- Unfortunately, this distinction is not theoretically consistent with 
the physical case m = rriyai = '^risea- 

In this work, I will describe an efficient method for simulating lattice QCD with dynam- 
ical fermions, which is particularly useful for investigating the thermodynamical properties 
and chiral properties. Results will be provided for the chiral condensate and gauge action 
in QCD on the 10'* lattice for flavor number Nj = 1, 2, 3, 4, and many values of quark mass 
m and coupling constant g. New results in the chiral limit for the gauge action and its 
correlation with the chiral condensate will also be presented. 



2 METHODOLOGY 

The method described here is a generahzation of the microcanonical fermionic average 
method [|T^] to QCD in 3+1 dimensions. Using the identity 



-^« = JdEdl^Y. ReTr(f/p) - 6^^] e*^^^^, 



e-*« = dE6\—Y^ ReTr(f/p) - 6VE e*^^^^, (4) 

one can rewrite the partition function as a one-dimensional integral 

Z= /rfEe-^^*(^'"^'^/'^). (5) 



Here S^^ is the full effective action as a function of plaquette energy E, quark mass m, 
flavor number Nf, and coupling f3 defined by 

S''^{E,m,Nf,(3) = -\nn{E)-Qf3VE 

+ Sf{E,m,Nf), (6) 

where V is the total number of lattice sites, 

n{E) = j[dU] 6 (^i- Y: ReTr(f/,) - 6VE^ (7) 

is the density of states at the given E, and 

Sf = -\n{{detA{m,U)ff^^)E 



''NcV/2 



Nf/4: 



M n (A?(f/)+rn^) )e (8) 



i=l 



is the effective fermionic action. Non-vanishing Sf^ implies the interaction between quarks 
and gluons. Here Aj(f/) is the i-th positive eigenvalue of the massless Dirac operator 
A(m = 0) and (...)_£; is the average of the observable over configurations with the probability 
distribution 6 fX]pR'6Tr(t/p)/A^c — QVEj /n{E) at fixed plaquette energy E. S'f does not 
depend on /?. Once I compute the the positive eigenvalues of A(m = 0) for pure gauge 
configurations at a set of E, I can obtain at no cost Sf^ and therefore the partition function 
for any m, Nf and /?. Then the thermodynamical properties can be obtained from the 
derivatives of the partition function. For example the chiral condensate and the vacuum 
expectation value of the plaquette are 

(ReTiUp) _ 1 d\nZ 



/ciEEe-^^"(^''"'^/'^) 



i^i;) 



Z 
dlnZ 



VN^Nf/A dm 



_1 / dE ^e-^^''(^-"^'^/'^) 

■^ rim. 



dm 



VN^Nf/A Z 



VNrZ 

(9) 



{T!:£(",^ri^{det/\{m^U)r^")E 



{{det^{m,U)f''^) 



E 



We can also calculate the specific heat, chiral susceptibility and other local quantities. One 
prominent feature is that the effective action and partition function are also calculable in 
the chiral limit, which is very useful for studying the chiral properties. 

3 IMPLEMENTATION 

The basic idea of the algorithm described in Sect. |^ is to compute the effective action in 
Eq. d^). The density of states n{E) can be directly evaluated by numerically integrating 
out the quenched SU(3) data: 



\\in{E) 



rE 

6 / dE' f3{E', Nf = 0) + const.. (10) 

Jo 



V 

cff 



The most important and time consuming part of the work, is the computation of Sf 
in Eq. (|^) by averaging out the fermionic determinant over configurations at the given 
E. These configurations can be generated by microcanonical or over-relaxation processes. 



For SU(3), although there exists a microcanonical algorithm |jT6|, it is very difficult to 
implement it ergodically. To solve the problem, I use a different prescription, over-relation- 
updating the subgroups SU(2) of SU(3) described as follows. Let U be an SU(3) link to be 
updated, and R the sum of 6 staples. 



(a) Find the 2x2 block of UR by striking out the 3rd row and column ||17|| : 




(11) 
\ y< X X J 

Write 

B = ro + zf.a=( ^' I'], (12) 



with ro and r being complex. Then find B' , and the norm k: 



B' 



Re(ro) + iRe(r) ■ a 

( {Br + Bt)/2 {B,-B*)/2 
[ {B, - Bl)/2 {B\ + B,)/2 



k = Vdet B' 
so that B' /k is a SU(2) matrix. Denote 

c = {{By/k) 

and 

oi = 



^^2 


/Ci 


C3 


"J - 


"U^ 


C^4 


Ci 


Cs 


\ 


C2 


C4 







1 


/ 



Perform over-relaxation updates 113 of the SU(2) subgroups: U' = aiU. 
(b) Find the 2x2 block of U'R by striking out the 2nd row and column: 




Repeat the procedure of finding C as Eqs. (|T2D? (HID, (|Tj) and denote 



a2 



Ci C3 
1 

C2 C4 



(c) Find the 2x2 block of U" R by striking out the 1st row and column: 




Repeat the procedure of finding C as Eqs. (|T2D? (|T3|), (115) and denote 



as 



1 
Ci C3 

C2 C4 



(13) 



(14) 



(15) 



(16) 



(17) 



:i8) 



(19) 



The new link is now Unew = asU" . The processes above satisfy Tr(f/„etu-R) = Tr(f/"_R) 
Tr{U'R) = Tt{UR), where the plaquette energy remains unchanged. 



The Lanczos algorithm [|19| was designed for calculating the eigenvalues of a large sparse 
matrix, but the rounding errors grow exponentially for the larger eigenvalues. I use the 
modified Lanczos algorithm ^^ ^ to solve this problem so that all the true eigenvalues 



can be found. On a d dimensional lattice with V lattice sites, the following sum rules 



NcV/2 

i=l 

NcV/2 

E Af = 



Ad -I , , 
-^^-{d-l)E 



-^- (20) 



can be used to check the accuracy of the eigenvalues of A(m = 0). 

4 RESULTS 

4.1 Microcanonical updates 

I describe now the QCD data on the 10^ lattice, obtained on a workstation. For the 
gauge fields, periodic boundary conditions are used. For fermions anti-periodic boundary 
condition in the time direction is implemented. 

First, I calculate the quenched mean plaquette energy as a function of /3, using the 
Cabibbo-Marinari |]T^ algorithm. High precision for this quantity can be easily achieved. 



Figure |] shows the result. Then I calculate the density of states using Eq. (0), by 
numerically integrating out the interpolated data for /3(-E) according to the trapezoidal 
rule. The result is shown in Fig. |], while the irrelevant constant in Eq. ([T0|) is ignored. 

At each fixed E, 10000-40000 pure gauge configurations are generated, each one is 
separated by 100 over-relation updates. 100-400 de-correlated configurations are used for 
the diagonalization of the massless Kogut-Susskind fermionic matrix A(m = 0). Their 
eigenvalues are stored on a disk with double precision, and the relative errors of the sum 
rules in Eq. (|20D are the order of 10~® and 10~^ respectively. Figure |^ plots the quenched 



chiral condensate 






versus m at P = 6.05. The results are in good agreement with Chen's data |2^, although 
Chen used the fermionic matrix inversion method on a much larger lattice 16^ x 32. This 
provides an additional check of the eigenvalues. 

From the eigenvalues of A{m = 0) at 16 values of -E G [0,1), we can compute the 
effective fermionic action in Eq. (j^) at any m and Nf at almost no cost. Figures §,^,P, and 
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1^ show the effective fermionic action (norinahzed by the volume) versus E for Nf = 1,2,3 
and 4 respectively. I calculated the effective fermionic action for 15 bare quark masses in 
m G [0, 0.1], but only illustrate 3 of them. Statistical errors are less than 0{1/V) = 0(10^^) 
and invisible at the scale of the figures. The shapes of the curves are quite similar, but the 
scales are quite different. The slope is small for small E < 0.4 and large E > 0.7, but it 
is big for intermediate E G (0.4,0.7). We find S'f oc Nf, which means the effects of the 
dynamical quarks are proportional to the flavor number, and are significant for intermediate 
E. 

With the density of states and the effective fermionic action, we can construct the full 
effective action in Eq. (^ as a function of E for any given Nf, m and (3, using the Newton 
polynomial interpolation. 

4.2 Themodynamical observables with sea quarks 

The thermodynamical quantities can be obtained by numerically integrating out the one- 
dimensional integrals in Eq. (|^). 

The mean plaquette energy Ep versus P is plotted in Fig. || for m = nisea = 0.01 and 
different Nf. For Nf = 2 and p = 5.7, Ep = 0.5771624 ± 6.2251091 x 10-^ consistent with 
Chen's HMD data Ep = 0.577386(17) on the 16^ x 40 in Ref. H. For Nf = A and /3 = 5.4, 



Ep = 0.5642520 ± 8.0579519 x 10-^ consistent with Chen's HMC data Ep = 0.560334(15) 
on the 16^ x 32 in Ref. [^]. This also means that finite size effects on Ep are not significant. 
As seen in this figure, at some fixed /3, Ep increases with Nf. The sea quark effects become 
most important around Ep ^ 0.5. We can understand this phenomenon by looking at Figs. 
^^,^ and ^ from which one observes maximum slope around E ^ 0.5. The lighter the 
quark mass, the bigger the slope is. Therefore, we have a mechanism: the sea quark effects 
reach their maximum where the effective fermionic action versus the gauge energy has a 
maximum slope. 

Figure ^ shows the data for Ep in the chiral limit m = rrisea = and different Nf. This 
quantity has been very useful for analyzing the chiral properties of the QED system ||23|| . 
We believe it will also be useful for QCD. 

Figure |10| depicts < ipip > versus m at f3 = 5.5 for different Nf. The data from the 
Langevin algorithm [^ on the 8^ x 18 lattice are also included for comparison. (The minor 
difference between our Nf = 2 data and those in |2^] might be accounted by different lattice 



sizes used or systematic error involved). Again, one sees < ipip > decreases significantly 
with Nf. 

4.3 Correlation between the chiral condensate and gauge action 
in the chiral limit 

The correlation function between the plaquette and the chiral condensate indicates the 
interaction between sea quarks and gluons. In the quenched case, it is zero, while for full 

7 



QCD, it is not. It can be computed by measuring the mass derivative of the plaquette 
mean value, since 

dm A \ N^ ^^^' 7 ^ ' 

where 

^^ = p4:;TrA-^ (23) 



A direct calculation of the r.h.s. of Eq. (22) needs very high statistics. Also it is impossible 
to calculate the r.h.s. when m = Q because on a finite lattice {ipip) vanishes identically in 
this limit. However, the l.h.s. of Eq. (^) can easily be calculated in our method. Figure 
|TT] plots the results Ep versus m for different j3 and Nf. The mass derivative is obtained 
by a linear fit to the data, which results at m = are given in Tab. I. 

5 OUTLOOK 



In this paper, I have extended the microcanonical fermionic average method [0 to QCD in 
3+1 dimensions, which allows us to search the parameter space {Nf, (3, m) with much lower 
computational cost. I have provided new data for the mean value of chiral condensates 
< 'ijjijj >, plaquette energy Ep, and their correlation function, including the results in the 
chiral limit for the latter two quantities. 

The disadvantages of the method is that it is not easy to calculate physical observables 
beyond the thermodynamical quantities (e.g. the spectrum). The storage of the eigenvalues 
of the fermionic matrix needs big hard disk space. Although the method is applicable to 
the chiral limit, the systematic errors become larger when the flavor number and /3 are 
larger, and the quark mass is smaller. 

As the algorithm |T5| has been useful for analyzing the QED phase structure, I believe 
the method developed in this work will also be useful for QCD. Since the algorithm also 



works in the chiral limit, the study of spontaneous chiral symmetry breaking [^ will be 
an interesting application of this work and will be reported elsewhere. 
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Table 1: Mass derivative of the mean plaquette energy at m = for /3 = 5.2. 



Nf 


dEp/dm 








1 


-0.0357841 ± 0.0003783 


2 


-0.0815237 ±0.001083 


3 


-0.327659 ±0.0182 


4 


-0.306373 ±0.000517 
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Figure 1: The vacuum expectation value of the plaquette energy from quenched simulation. 
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Figure 2: — ln[n(-E)]/y as a function of E. 
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Figure 3: The quenched chiral condensate < ipip > versus m ai f3 = 6.05. 
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Figure 4: —Sj/V as a function of i? for Nf = 1. 
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Figure 5: —Sf/V as a function of E for Nj = 2. 
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Figure 6: —Sj/V as a function of E for Nf = 3. 
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Figure 7: —Sf/V as a function of i? for Nf = 4. 
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Figure 8: The mean plaquette energy Ep versus (3 at m = 0.01. 
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Figure 9: The mean plaquette energy Ep versus P at m = 0, obtained by the method 
described in this paper. 
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Figure 10: Chiral condensate < ipip > versus m at /5 = 5.5. 
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Figure 11: Ep versus m aX j3 = 5.2. 
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